Spatial data are used across a wide range of fields to support decision-making, including environment, public health, ecology, agriculture, urban planning, economy, and society.
Areal data usually arise when the number of events corresponding to some variable of interest are aggregated in areas.
library(sf)Warning: package 'sf' was built under R version 4.4.3Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUElibrary(mapview)
d <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
mapview(d, zcol = "SID74")Example of areal data. Number of sudden infant deaths in counties of North Carolina, USA, in 1974.
library(spData)Warning: package 'spData' was built under R version 4.4.2library(ggplot2)
d <- st_read(system.file("shapes/columbus.gpkg", package = "spData"), quiet = TRUE)
st_crs(d) <- NA
ggplot(d) + geom_sf(aes(fill = INC))
Example of areal data. Household income in $1000 USD in neighborhoods in Columbus, Ohio, in 1980.
library(terra)Warning: package 'terra' was built under R version 4.4.3terra 1.8.60d <- rast(system.file("ex/elev.tif", package = "terra"))
plot(d)
Example of areal data. Elevation at raster grid cells covering Luxembourg.
For example, we can use air pollution measurements at a number of monitoring stations to predict air pollution at other locations taking into account spatial autocorrelation and other factors that are known to predict the outcome of interest.
library(sp)
library(sf)
library(mapview)
data(meuse)
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
mapview(meuse, zcol = "lead", map.types = "CartoDB.Voyager")Example of geostatistical data. Topsoil lead concentrations at locations sampled in a flood plain of the river Meuse, The Netherlands.
library(spData)
mapview(properties, zcol = "prpsqm")Example of geostatistical data. Price per square meter of a set of apartments in Athens, Greece, in 2017.
library(malariaAtlas)Warning: package 'malariaAtlas' was built under R version 4.4.3Registered S3 method overwritten by 'malariaAtlas':
method from
autoplot.default ggplot2d <- getPR(country = "Zimbabwe", species = "BOTH")Loading ISO 19139 XML schemas...Loading ISO 19115-3 XML schemas...Loading ISO 19139 codelists...Please Note: Because you did not provide a version, by default the version being used is 202406 (This is the most recent version of PR data. To see other version options use function listPRPointVersions)Warning: The `sourcedata` argument of `isAvailable_pr()` is deprecated as of
malariaAtlas 1.6.0.
ℹ The argument 'sourcedata' has been deprecated. It will be removed in the next
version. It has no meaning.
ℹ The deprecated feature was likely used in the malariaAtlas package.
Please report the issue at
<https://github.com/malaria-atlas-project/malariaAtlas/issues>.Please Note: Because you did not provide a version, by default the version being used is 202406 (This is the most recent version of PR data. To see other version options use function listPRPointVersions)Confirming availability of PR data for: Zimbabwe...PR points are available for Zimbabwe.Attempting to download PR point data for Zimbabwe ...Data downloaded for Zimbabwe.NOTE: All available data for this query was downloaded for both species,
but there are no PR points for P. vivax in this region in the MAP database.
To check endemicity patterns or to contribute data, visit malariaatlas.org OR email us at map@bdi.ox.ac.uk.ggplot2::autoplot(d)Warning: The `format` argument of `getShp()` is deprecated as of malariaAtlas 1.6.0.
ℹ The argument 'format' has been deprecated. It will be removed in the next
version. Admin boundaries will be correctly plotted using autoplot without
the argument.
ℹ The deprecated feature was likely used in the malariaAtlas package.
Please report the issue at
<https://github.com/malaria-atlas-project/malariaAtlas/issues>.Please Note: Because you did not provide a version, by default the version being used is 202403 (This is the most recent version of admin unit shape data. To see other version options use function listShpVersions)
Example of geostatistical data. Malaria prevalence at specific locations in Zimbabwe.
Point patterns arise when the variable to be analyzed corresponds to the location of events.
library(spatstat)Warning: package 'spatstat' was built under R version 4.4.2Loading required package: spatstat.dataWarning: package 'spatstat.data' was built under R version 4.4.2Loading required package: spatstat.univarWarning: package 'spatstat.univar' was built under R version 4.4.2spatstat.univar 3.1-1Loading required package: spatstat.geomWarning: package 'spatstat.geom' was built under R version 4.4.2spatstat.geom 3.3-5
Attaching package: 'spatstat.geom'The following objects are masked from 'package:terra':
area, delaunay, is.empty, rescale, rotate, shift, where.max,
where.minLoading required package: spatstat.randomWarning: package 'spatstat.random' was built under R version 4.4.2spatstat.random 3.3-2Loading required package: spatstat.exploreWarning: package 'spatstat.explore' was built under R version 4.4.2Loading required package: nlmespatstat.explore 3.3-4Loading required package: spatstat.modelWarning: package 'spatstat.model' was built under R version 4.4.2Loading required package: rpartspatstat.model 3.3-4Loading required package: spatstat.linnetWarning: package 'spatstat.linnet' was built under R version 4.4.2spatstat.linnet 3.2-5
spatstat 3.3-1
For an introduction to spatstat, type 'beginner' plot(clmfires, use.marks = FALSE, pch = ".")
Data clmfires is a marked point pattern containing information of each fire, depicts the location of the fires without the mark.
library(spatstat)
plot(hamster)
This figure also shows the positions of cell nuclei in a histological section of a tissue from a lymphoma in the kidney of a hamster from spatstat. The nuclei are classified as either “pyknotic” (corresponding to dying cells) or “dividing” (corresponding to cells arrested in the act of dividing).
library(sparr)Warning: package 'sparr' was built under R version 4.4.3
Welcome to
_____ ___ ____ ____ ____
/ ___// _ \/ _ \/ __ \/ __ \
\__ \/ ___/ __ / ___/ ___/
___/ / / / / / / /\ \/ /\ \
/____/_/ /_/ /_/_/ \__/ \_\ v2.3-16
- type news(package="sparr") for an overview
- type help("sparr") for documentation
- type citation("sparr") for how to citedata(pbc)
plot(unmark(pbc[which(pbc$marks == "case"), ]), main = "cases")
axis(1)
axis(2)
title(xlab = "Easting", ylab = "Northing")
plot(unmark(pbc[which(pbc$marks == "control"), ]), pch = 3, main = "controls")
axis(1)
axis(2)
title(xlab = "Easting", ylab = "Northing")
Examples of point patterns. Top: Locations of fires in Castilla-La Mancha, Spain, between 1998 and 2007. Bottom: Locations and types of cells in a tissue.
Spatio-temporal data arise when information is both spatially and temporally referenced. Thus, we can consider spatial data as temporal aggregations or temporal snapshots of a spatio-temporal process.
Spatial functional data arise when the three types of spatial data (areal, geostatistical, and point patterns) are combined with random functions.
Besides the three classical types of spatial data (i.e., areal, geostatistical, and point patterns), we can also consider other spatial data such as flows containing the number of individuals or other elements moving between locations.
library("epiflows")Warning: package 'epiflows' was built under R version 4.4.3epiflows is loaded with the following global variables in `global_vars()`:
coordinates, pop_size, duration_stay, first_date, last_date, num_casesdata("Brazil_epiflows")
loc <- merge(x = YF_locations, y = YF_coordinates,
by.x = "location_code", by.y = "id", sort = FALSE)
ef <- make_epiflows(flows = YF_flows, locations = loc,
coordinates = c("lon", "lat"),
pop_size = "location_population",
duration_stay = "length_of_stay",
num_cases = "num_cases_time_window",
first_date = "first_date_cases",
last_date = "last_date_cases")
ef
/// Epidemiological Flows //
// class: epiflows, epicontacts
// 15 locations; 100 flows; directed
// optional variables: coordinates, pop_size, duration_stay, num_cases, first_date, last_date
// locationsWarning: `as_data_frame()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` (with slightly different semantics) to convert to a
tibble, or `as.data.frame()` to convert to a data frame.
ℹ The deprecated feature was likely used in the epiflows package.
Please report the issue at <https://github.com/reconhub/epiflows/issues>.# A tibble: 15 × 8
id location_population num_cases_time_window first_date_cases
<chr> <dbl> <dbl> <fct>
1 Espirito Santo 3973697 2600 2017-01-04
2 Minas Gerais 20997560 4870 2016-12-19
3 Rio de Janeiro 16635996 170 2017-02-19
4 Sao Paulo 44749699 200 2016-12-17
5 Southeast Brazil 86356952 7840 2016-12-17
6 Argentina NA NA <NA>
7 Chile NA NA <NA>
8 Germany NA NA <NA>
9 Italy NA NA <NA>
10 Paraguay NA NA <NA>
11 Portugal NA NA <NA>
12 Spain NA NA <NA>
13 United Kingdom NA NA <NA>
14 United States of … NA NA <NA>
15 Uruguay NA NA <NA>
# ℹ 4 more variables: last_date_cases <fct>, length_of_stay <dbl>, lon <dbl>,
# lat <dbl>
// flows
# A tibble: 100 × 3
from to n
<chr> <chr> <dbl>
1 Espirito Santo Italy 2828.
2 Minas Gerais Italy 15714.
3 Rio de Janeiro Italy 8164.
4 Sao Paulo Italy 34039.
5 Southeast Brazil Italy 76282.
6 Espirito Santo Spain 3270.
7 Minas Gerais Spain 18176.
8 Rio de Janeiro Spain 9443.
9 Sao Paulo Spain 39371.
10 Southeast Brazil Spain 88231.
# ℹ 90 more rowsWe can use this data to create an epiflows object called ef that allows us to use the prediction and visualization functions.
The sf package allows us to work with vector data which is used to represent points, lines, and polygons.
Vector data are commonly stored in the shapefile format.
A shapefile is not a single file, but a collection of related files that work together.
The three mandatory components are:
.shp – contains the geometry data,
.shx – a positional index that allows quick access to the geometry,
.dbf – stores attribute data for each feature.
Additional optional files may include:
.prj – specifies the map projection,
.sbn and .sbx – spatial index files,
.shp.xml – contains spatial metadata in XML format.
library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")map<-st_read(pathshp,quiet= TRUE)
class(map)[1] "sf" "data.frame"head(map)Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS: NAD27
AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
NWBIR74 BIR79 SID79 NWBIR79 geometry
1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...
6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3...plot(map[1]) #plotfirstattribute
Map of the first attribute of the ’sf‘ object representing the counties of North Carolina, USA.
Raster data (also referred to as grid data) is a spatial data structure that divides the region of study into rectangles of the same size called cells or pixels, and that can store one or more values for each of these cells.
Raster data often come in GeoTIFF format which has extension .tif.
library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- terra::rast(pathraster)
rclass : SpatRaster
size : 90, 95, 1 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : elev.tif
name : elevation
min value : 141
max value : 547 plot(r)
Map of raster data representing the elevation of Luxembourg.
In a geographic CRS, latitude and longitude values are used to identify locations on the Earth’s three-dimensional ellipsoid surface.
Latitude values range from –90° at the South Pole to +90° at the North Pole. Longitude values measure the angle east or west of the Prime Meridian.
Latitude and longitude coordinates may be expressed in degrees, minutes, and seconds, or in decimal degrees.
Projected CRSs use Cartesian coordinates to reference a location on a two-dimensional representation of the Earth.
A common projection is the Universal Transverse Mercator (UTM) projection. This projection is conformal, meaning it preserves angles and therefore maintains shapes over small areas.
Most common CRSs can be specified using their EPSG (European Petroleum Survey Group) codes or their Proj4 strings. Common spatial projections can be found at https://spatialreference.org/ref/. Details of a given projection can be inspected using the st_crs() function of the sf package. For example, EPSG code 4326 refers to the WGS84 longitude/latitude projection.
st_crs("EPSG:4326")$Name[1] "WGS 84" st_crs("EPSG:4326")$proj4string[1] "+proj=longlat +datum=WGS84 +no_defs"st_crs("EPSG:4326")$epsg[1] 4326library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
map <- st_read(pathshp, quiet = TRUE)
# Get CRS
# st_crs(map)
# Transform CRS
map2 <- st_transform(map, crs = "EPSG:4326")
# Get CRS
# st_crs(map2)library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
# Get CRS
# crs(r)
# Transform CRS
r2 <- terra::project(r, "EPSG:2169")
# Get CRS
# crs(r2)Alternatively, if we want the transformed data to align exactly with other raster data we are using, we can project it using an existing raster that has the desired geometry.
# x is existing raster
# r is raster we project
# r2 <- terra::project(r, x)Before the sf package was developed, the sp package was used to represent and work with vector spatial data. The sp, rgdal (Bivand et al., 2023), rgeos (Bivand and Rundel, 2022), and maptools (Bivand and Lewin-Koh, 2022) packages are no longer maintained and have been retired.
library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(pathshp, quiet = TRUE)
class(nc)[1] "sf" "data.frame"A sf object contains the following objects of class sf, sfc and sfg: • sf (simple feature): each row of the data.frame is a single simple feature consisting of attributes and geometry. • sfc (simple feature geometry list-column): the geometry column of the data.frame is a list-column of class sfc with the geometry of each simple feature. • sfg (simple feature geometry): each of the rows of the sfc list-column corresponds to the simple feature geometry (sfg) of a single simple feature.
head(nc)Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS: NAD27
AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
NWBIR74 BIR79 SID79 NWBIR79 geometry
1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...
6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3...plot(nc)Warning: plotting the first 10 out of 14 attributes; use max.plot = 14 to plot
all
object representing the counties of North Carolina, USA, and associated information.
nc[1,] # first rowSimple feature collection with 1 feature and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
Geodetic CRS: NAD27
AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1 10
BIR79 SID79 NWBIR79 geometry
1 1364 0 19 MULTIPOLYGON (((-81.47276 3...nc[nc$NAME == "Ashe", ] # row with NAME "Ashe"Simple feature collection with 1 feature and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
Geodetic CRS: NAD27
AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1 10
BIR79 SID79 NWBIR79 geometry
1 1364 0 19 MULTIPOLYGON (((-81.47276 3...nc[1, "NWBIR74"] # first row, column with name NWBIR74Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
Geodetic CRS: NAD27
NWBIR74 geometry
1 10 MULTIPOLYGON (((-81.47276 3...nc[1, "NWBIR74", drop = TRUE] # drop geometry[1] 10
attr(,"class")
[1] "numeric"# Geometries printed in abbreviated form
st_geometry(nc)Geometry set for 100 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS: NAD27
First 5 geometries:MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...# View complete geometry by selecting one
st_geometry(nc)[[1]]MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436)))Simple feature geometries (sfg objects) can be, for example, of type POINT (single point), MULTIPOINT (set of points), or POLYGON (polygon), and can be created using the functions st_point(), st_multipoint(), and st_polygon(), respectively. Here, we create an sf object containing two single points, a set of points, and a polygon, each with one attribute.
library(sf)
library(ggplot2)
# Single points (point as a vector)
p1_sfg <- st_point(c(2, 2))
p2_sfg <- st_point(c(2.5, 3))
# Set of points (points as a matrix)
p <- rbind(
c(6, 2),
c(6.1, 2.6),
c(6.8, 2.5),
c(6.2, 1.5),
c(6.8, 1.8)
)
mp_sfg <- st_multipoint(p)
# Polygon: sequence of points forming a closed, non-self-intersecting ring
# The first ring denotes the exterior ring;
# zero or more subsequent rings denote holes within it
p1 <- rbind(c(10, 0), c(11, 0), c(13, 2),
c(12, 4), c(11, 4), c(10, 0))
p2 <- rbind(c(11, 1), c(11, 2), c(12, 2), c(11, 1))
pol_sfg <- st_polygon(list(p1, p2))
# Create sf object
p_sfc <- st_sfc(p1_sfg, p2_sfg, mp_sfg, pol_sfg)
df <- data.frame(v1 = c("A", "B", "C", "D"))
p_sf <- st_sf(df, geometry = p_sfc)
# Plot single points, multipoint, and polygon
ggplot(p_sf) +
geom_sf(aes(col = v1), size = 3) +
theme_bw()
sf object representing two single points, a set of points, and a polygon, with one attribute.
Common functions to manipulate sf objects include:
• st_read() – reads an sf object
• st_write() – writes an sf object
• st_crs() – gets or sets a new coordinate reference system (CRS)
• st_transform() – transforms data to a new CRS
• st_intersection() – intersects sf objects
• st_union() – combines several sf objects into one
• st_simplify() – simplifies an sf object
• st_coordinates() – retrieves coordinates of an sf object
• st_as_sf() – converts another object type to an sf object
We can read a sf object as follows:
library(sf)
library(ggplot2)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))We can inspect the first rows of the sf object map with head().
head(map)Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS: NAD27
# A tibble: 6 × 15
AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
<dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <dbl> <dbl>
1 0.114 1.44 1825 1825 Ashe 37009 37009 5 1091 1 10
2 0.061 1.23 1827 1827 Alleg… 37005 37005 3 487 0 10
3 0.143 1.63 1828 1828 Surry 37171 37171 86 3188 5 208
4 0.07 2.97 1831 1831 Curri… 37053 37053 27 508 1 123
5 0.153 2.21 1832 1832 North… 37131 37131 66 1421 9 1066
6 0.097 1.67 1833 1833 Hertf… 37091 37091 46 1452 7 954
# ℹ 4 more variables: BIR79 <dbl>, SID79 <dbl>, NWBIR79 <dbl>,
# geometry <MULTIPOLYGON [°]>We can delete some of the polygons by taking a subset of the rows of map.
# Delete polygon
map <- map[-which(map$FIPS %in% c("37125", "37051")), ]
ggplot(map) + geom_sf(aes(fill = SID79))
We can use st_union() with argument by_feature = FALSE to combine all geometries together.
# Combine geometries
ggplot(st_union(map, by_feature = FALSE) %>%
st_sf()) + geom_sf()although coordinates are longitude/latitude, st_union assumes that they are
planar
The boundaries of a map can be simplified with the st_simplify() function.
# Simplify
map_utm <- st_transform(map, 32617)
map_simplified <- st_simplify(map_utm, dTolerance = 10000)
ggplot(map_simplified) +
geom_sf(fill = "lightblue", color = "gray40") +
theme_minimal() +
labs(title = "Simplified North Carolina Map (10 km tolerance)")
The st_as_sf() function allows us to convert a foreign object into an sf object. For example, we can have a data frame containing the coordinates of several locations along with their attributes, which we can then convert into an sf object.
library(sf)
library(mapview)
d <- data.frame(
place = c("London", "Paris", "Madrid", "Rome"),
long = c(-0.118092, 2.349014, -3.703339, 12.496366),
lat = c(51.509865, 48.864716, 40.416729, 41.902782),
value = c(200, 300, 400, 500))
class(d)[1] "data.frame"dsf <- st_as_sf(d, coords = c("long", "lat"))
st_crs(dsf) <- 4326
class(dsf)[1] "sf" "data.frame"mapview(dsf)We can use the st_intersects() function from sf to count the number of points within the polygons of an sf object. The returned object is a list containing the feature IDs that intersect each polygon. We can then use the lengths() function to calculate the number of points inside each feature.
library(sf)
library(ggplot2)
# Map with divisions (sf object)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))
# Generate random points within the polygons (sf geometry list-column sfc)
points <- st_sample(map, size = 100, exact = FALSE)although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar# Convert points (sfc) into an sf object so it can be used in ggplot easily
points_sf <- st_sf(geometry = points)
# Plot map and points
ggplot()+geom_sf(data= map)+geom_sf(data=points)
Points within polygons.
# Intersection: which points fall in which polygons
inter <- st_intersects(map, points_sf)although coordinates are longitude/latitude, st_intersects assumes that they
are planar# Add point count to each polygon
map$count <- lengths(inter)
# Plot number of points within each polygon
ggplot(map)+geom_sf(aes(fill = count))
Number of points within polygons.
library(sf)
library(ggplot2)
library(dplyr)
Attaching package: 'dplyr'The following object is masked from 'package:nlme':
collapseThe following objects are masked from 'package:terra':
intersect, unionThe following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, union# Map with divisions (sf object)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))
# Generate random points within the polygons
points <- st_sample(map, size = 3) %>%
st_as_sf()although coordinates are longitude/latitude, st_intersects assumes that they
are planaralthough coordinates are longitude/latitude, st_intersects assumes that they
are planar# Find which polygon each point falls inside
# (note: first argument = points, second = polygons)
inter <- st_intersects(points, map)although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar# Add column with the name of the area containing each point
points$areaname <- map[unlist(inter), "NAME", drop = TRUE]
# Show points with corresponding area name
pointsSimple feature collection with 3 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -79.0806 ymin: 34.94656 xmax: -76.03757 ymax: 35.96621
Geodetic CRS: NAD27
x areaname
1 POINT (-79.0806 34.94656) Hoke
2 POINT (-76.03757 35.91925) Tyrrell
3 POINT (-78.32428 35.96621) Franklin# Map
ggplot(map) +
geom_sf() +
geom_sf(data = points) +
geom_sf_label(
data = map[unlist(inter), ],
aes(label = NAME),
nudge_y = 0.2
)Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
Sometimes, a map and its corresponding data are available separately, and we may want to create an sf object that combines the map with its related data for easier manipulation and plotting. This can be done by joining the map and data using the left_join() function from the dplyr package.
library(rnaturalearth)Warning: package 'rnaturalearth' was built under R version 4.4.3map <- ne_countries(returnclass = "sf")library(wbstats)Warning: package 'wbstats' was built under R version 4.4.3indicators <- wb_search(pattern = "pollution")
d <- wb_data(indicator = "EN.ATM.PM25.MC.M3", start_date = 2016, end_date = 2016)Then, we use the left_join() function from dplyr to merge the map and data, specifying the variable(s) to join by using the by argument. In this example, we use the ISO3 standard country codes instead of country names, as names may differ between the map and the data frame. Figure 3.7 shows the resulting map created using ggplot2.
library(dplyr)
library(ggplot2)
library(viridis)Loading required package: viridisLitemap1 <- left_join(map, d, by = c("iso_a3" = "iso3c"))
ggplot(map1) + geom_sf(aes(fill = EN.ATM.PM25.MC.M3)) +
scale_fill_viridis() + labs(fill = "PM2.5") + theme_bw()
PM2.5 values in each of the world countries in 2016.
Note that when using left_join(), the resulting object’s class matches that of the first argument.
Thus, left_join(sf_object, data.frame_object) produces an sf object, while left_join(data.frame_object, sf_object) produces a data frame.
map1 <- left_join(map, d, by = c("iso_a3" = "iso3c"))
class(map1)[1] "sf" "data.frame"d1 <- left_join(d, map, by = c("iso3c" = "iso_a3"))
class(d1)[1] "tbl_df" "tbl" "data.frame"library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)
We can also use rast() to create a SpatRaster object by specifying the number of columns, number of rows, and the minimum and maximum x and y values.
r <- rast(ncol = 10, nrow = 10, xmin = -150, xmax = -80, ymin = 20, ymax = 60)
rclass : SpatRaster
size : 10, 10, 1 (nrow, ncol, nlyr)
resolution : 7, 4 (x, y)
extent : -150, -80, 20, 60 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) nrow(r) # number of rows[1] 10ncol(r) # number of columns[1] 10dim(r) # dimension[1] 10 10 1ncell(r) # number of cells[1] 100values(r) <- 1:ncell(r)
r2 <- r * r
s <- c(r, r2)
plot(s[[2]]) # layer 2
Many generic functions can be used to operate with rasters, as shown in the following examples.
plot(min(s))
plot(r + r + 10)
plot(round(r))
plot(r == 1)
pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)# Longitude and latitude values
long <- c(-0.118092, 2.349014, -3.703339, 12.496366)
lat <- c(51.509865, 48.864716, 40.416729, 41.902782)
longlat <- cbind(long, lat)
# CRS
crspoints <- "+proj=longlat +datum=WGS84"
# Attributes for points
d <- data.frame(
place = c("London", "Paris", "Madrid", "Rome"),
value = c(200, 300, 400, 500)
)
# SpatVector object
pts <- vect(longlat, atts = d, crs = crspoints)
pts class : SpatVector
geometry : points
dimensions : 4, 2 (geometries, attributes)
extent : -3.703339, 12.49637, 40.41673, 51.50986 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs
names : place value
type : <chr> <num>
values : London 200
Paris 300
Madrid 400# Plot
plot(pts)
library(terra)
r <- geodata::worldclim_country(country = "Spain", var = "tavg", res = 10, path = tempdir())
plot(r)
We can average the temperature raster data over months with the mean() function.
r <- mean(r)
plot(r)
We also download the map of Spain with the rnaturalearth package (South, 2017), and delete the Canary Islands region.
# Map
library(ggplot2)
map <- rnaturalearth::ne_states("Spain", returnclass = "sf")
map <- map[-which(map$region == "Canary Is."), ] # delete region
ggplot(map) + geom_sf()
We obtain the spatial extent of the map with terra:ext(). Then, we can use crop() to remove the part of the raster that is outside the spatial extent.
# Cropping
sextent <- terra::ext(map)
r <- terra::crop(r, sextent)
plot(r)
We can use the mask() function to convert all values outside the map to NA.
# Masking
r <- terra::mask(r, vect(map))
plot(r)
The aggregate() function in terra can be used to group raster cells and create a new raster with lower resolution (i.e., larger cells).
# Aggregating
r <- terra::aggregate(r, fact = 20, fun = "mean", na.rm = TRUE)
plot(r)
library(terra)
# Raster (SpatRaster)
r <- rast(system.file("ex/elev.tif", package = "terra"))
# Polygons (SpatVector)
v <- vect(system.file("ex/lux.shp", package = "terra"))We use the terra functions centroids() to get the centroids of the division polygons, and crds() to obtain their coordinates.
points <- crds(centroids(v))plot(r)
plot(v, add = TRUE)
points(points)
Elevation raster, and division and centroids of polygons in Luxembourg.
We can obtain raster values at points using extract(), with the first argument as the SpatRaster object and the second argument as the data frame containing the points.
# data frame with the coordinates
points <- as.data.frame(points)
valuesatpoints <- extract(r, points)
cbind(points, valuesatpoints) x y ID elevation
1 6.009082 50.07064 1 444
2 6.127425 49.86614 2 295
3 5.886502 49.80014 3 382
4 6.165081 49.92886 4 404
5 5.914545 49.93892 5 414
6 6.378449 49.78511 6 320
7 6.311601 49.54569 7 193
8 6.346395 49.68742 8 228
9 5.963503 49.64159 9 313
10 6.023816 49.52331 10 282
11 6.167624 49.61815 11 328
12 6.113598 49.75744 12 221We can also use extract() to obtain raster values from SpatRaster objects within polygons of class SpatVector.
# Extracted raster cells within each polygon
head(extract(r, v, na.rm = TRUE)) ID elevation
1 1 547
2 1 485
3 1 497
4 1 515
5 1 515
6 1 515# Extracted raster cells and percentage of area covered within each polygon
head(extract(r, v, na.rm = TRUE, weights = TRUE)) ID elevation weight
1 1 NA 0.04545454
2 1 NA 0.10909091
3 1 529 0.24545454
4 1 542 0.46363635
5 1 547 0.68181816
6 1 535 0.11818181# Average raster values by polygon
v$avg <- extract(r, v, mean, na.rm = TRUE)$elevation
# Area-weighted average raster values by polygon (weights = TRUE)
v$weightedavg <- extract(r, v, mean, na.rm = TRUE, weights = TRUE)$elevationlibrary(ggplot2)
library(tidyterra)Warning: package 'tidyterra' was built under R version 4.4.3Registered S3 method overwritten by 'tidyterra':
method from
autoplot.SpatRaster malariaAtlas
Attaching package: 'tidyterra'The following object is masked from 'package:stats':
filter# Plot average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = avg)) + scale_fill_terrain_c()
# Plot area-weighted average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = weightedavg)) + scale_fill_terrain_c()
Average and area-weighted average of elevation values in each of the divisions of Luxembourg.
Maps allow us to effectively communicate spatial information. Here, we demonstrate how to create both static and interactive maps using several mapping packages, including ggplot2 (Wickham et al., 2022a), leaflet (Cheng et al., 2022a), mapview (Appelhans et al., 2022), and tmap (Tennekes, 2022).
library(sf)
nameshp <- system.file("shape/nc.shp", package = "sf")
d <- st_read(nameshp, quiet = TRUE)
d$vble <- d$SID74
d$vble2 <- d$SID79library(ggplot2)
library(viridis)
ggplot(d) + geom_sf(aes(fill = vble)) +
scale_fill_viridis() + theme_bw()
Alternatively, we can specify a device driver (e.g., png, pdf), print the plot, and then close the device using dev.off().
#png("plot.png")
#ggplot(d) + geom_sf(aes(fill = vble)) + scale_fill_viridis() + theme_bw()
#dev.off()library(plotly)
Attaching package: 'plotly'The following object is masked from 'package:ggplot2':
last_plotThe following object is masked from 'package:stats':
filterThe following object is masked from 'package:graphics':
layoutg <- ggplot(d) + geom_sf(aes(fill = vble))
ggplotly(g)The sf object passed to leaflet() must have a geographic coordinate reference system (CRS) using latitude and longitude (EPSG:4326). Here, we use the st_transform() function from sf to convert the data d, which has CRS EPSG:4267, to CRS EPSG:4326.
st_crs(d)$epsg[1] 4267d <- st_transform(d, 4326)library(leaflet)Warning: package 'leaflet' was built under R version 4.4.2pal <- colorNumeric(palette = "YlOrRd", domain = d$vble)
l <- leaflet(d) %>% addTiles() %>%
addPolygons(color = "white", fillColor = ~ pal(vble),
fillOpacity = 0.8) %>%
addLegend(pal = pal, values = ~vble, opacity = 0.8)
lWe can also use the addMiniMap() function to include an inset map.
l %>% addMiniMap()For example, to save the leaflet map l as .png we can proceed as follows:
# webshot::install_phantomjs()
# Saves map.html
# library(htmlwidgets)
# saveWidget(widget = l, file = "map.html")
# Takes a screenshot of the map.html created and saves it as map.png
# library(webshot)
# webshot::install_phantomjs()
# webshot(url = "map.html", file = "map.png")library(mapview)
mapview(d, zcol = "vble")Maps created with mapview can also be customized to include elements such as legends and background maps.
library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
mapview(d, zcol = "vble", map.types = "CartoDB.DarkMatter",
col.regions = pal, layer.name = "SDI")An inset map can also be added by using the addMiniMap() function of leaflet.
map1 <- mapview(d, zcol = "vble")
leaflet::addMiniMap(map1@map)library(leaflet.extras2)Warning: package 'leaflet.extras2' was built under R version 4.4.3library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
# common legend
at <- seq(min(c(d$vble, d$vble2)), max(c(d$vble, d$vble2)), length.out = 8)
m1 <- mapview(d, zcol = "vble", map.types = "CartoDB.Positron",
col.regions = pal, at = at)
m2 <- mapview(d, zcol = "vble2", map.types = "CartoDB.Positron",
col.regions = pal, at = at)
m1 | m2library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
# common legend
at <- seq(min(c(d$vble, d$vble2)), max(c(d$vble, d$vble2)), length.out = 8)
m1 <- mapview(d, zcol = "vble", map.types = "CartoDB.Positron", col.regions = pal, at = at)
m2 <- mapview(d, zcol = "vble2", map.types = "CartoDB.Positron", col.regions = pal, at = at)
m <- leafsync::sync(m1, m2)
mThe tmap package (Tennekes, 2022) allows us to create static and interactive maps composed of multiple shapes and layers, with various styles.
library(tmap)Warning: package 'tmap' was built under R version 4.4.2tmap_mode("plot")ℹ tmap mode set to "plot".tm_shape(d) + tm_polygons("vble")
library(maps)Warning: package 'maps' was built under R version 4.4.2
Attaching package: 'maps'The following object is masked from 'package:viridis':
unempd <- world.cities
# Select South Africa
d <- d[which(d$country.etc == "South Africa"), ]
# Transform data to sf object
d <- st_as_sf(d, coords = c("long", "lat"))
# Assign CRS
st_crs(d) <- 4326d$vble <- d$pop
d$size <- sqrt(d$vble)/100
ggplot(d) + geom_sf(aes(col = vble, size = size)) +
scale_color_viridis()
A leaflet map can be created using addCircles() specifying the radius and color of the points.
pal <- colorNumeric(palette = "viridis", domain = d$vble)
leaflet(d) %>%
addTiles() %>%
addCircles(
lng = st_coordinates(d)[, 1],
lat = st_coordinates(d)[, 2],
radius = ~sqrt(vble) * 10,
color = ~pal(vble),
popup = ~name
) %>%
addLegend(
pal = pal,
values = ~vble,
position = "bottomright"
)To create a map with mapview, we set the size of the points with cex = “size”.
d$size <- sqrt(d$vble)
mapview(d, zcol = "vble", cex = "size")Finally, we use tm_dots() to create a map with tmap.
tmap_mode("view")ℹ tmap mode set to "view".tm_shape(d) + tm_dots("vble", scale = sqrt(d$vble)/500, palette = "viridis")
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_dots()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).[v3->v4] `tm_tm_dots()`: migrate the argument(s) related to the scale of the
visual variable `size` namely 'scale' (rename to 'values.scale') to size.scale
= tm_scale_continuous(<HERE>).
ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
a list: 'size.scale = list(<scale1>, <scale2>, ...)'library(terra)
filename <- system.file("ex/elev.tif", package = "terra")
r <- rast(filename)
# Transform data to sf object
d <- st_as_sf(as.data.frame(r, xy = TRUE), coords = c("x", "y"))
# Assign CRS
st_crs(d) <- 4326
# Plot
ggplot(d) + geom_sf() +
geom_raster(data = as.data.frame(r, xy = TRUE),
aes(x = x, y = y, fill = elevation))
To use the leaflet and mapview packages, we transform the data from class terra to RasterLayer with the raster::brick() function.
library(raster)
Attaching package: 'raster'The following object is masked from 'package:plotly':
selectThe following object is masked from 'package:tidyterra':
selectThe following object is masked from 'package:dplyr':
selectThe following object is masked from 'package:nlme':
getDatarb <- raster::brick(r)
pal <- colorNumeric("YlOrRd", values(r), na.color = "transparent")
leaflet() %>%
addTiles() %>%
addRasterImage(rb, colors = pal, opacity = 0.8) %>%
addLegend(
pal = pal,
values = values(r),
title = "Elevation"
)We create a map with mapview setting the title with the argument layer as follows:
mapview(rb, layer = "elevation")We create a map with tmap by using tm_raster().
tm_shape(r) + tm_raster(title = "elevation", palette = "YlOrRd")── tmap v3 code detected ───────────────────────────────────────────────────────[v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
visual variable `col` namely 'palette' (rename to 'values') to col.scale =
tm_scale(<HERE>).
[v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
"brewer.yl_or_rd"#devtools::install_github("FlowmapBlue/flowmapblue.R")
library(flowmapblue)
locations <- data.frame(
id = c(1, 2, 3),
name = c("New York", "London", "Rio de Janeiro"),
lat = c(40.713543, 51.507425,-22.906241),
lon = c(-74.011219,-0.127738,-43.180244)
)The data frame flows has the number of people moving between origin and destination locations.
flows <- data.frame(
origin = c(1, 2, 3, 2, 1, 3),
dest = c(2, 1, 1, 3, 3, 2),
count = c(42, 51, 50, 40, 22, 42)
)Finally, we call the flowmapblue() function, passing the locations and flows data frames, along with the Mapbox access token, and specifying options such as clustering and animation. If mapboxAccessToken is set to NULL, the flows between locations will still be visualized, but without a background map.
#flowmapblue(locations, flows, mapboxAccessToken = NULL, clustering = TRUE, darkMode = TRUE, animation = FALSE)# devtools::install_github("ropensci/rnaturalearthhires")
library(rnaturalearth)
library(sf)
library(ggplot2)
library(viridis)
library(patchwork)
Attaching package: 'patchwork'The following object is masked from 'package:raster':
areaThe following object is masked from 'package:spatstat.geom':
areaThe following object is masked from 'package:terra':
areamap1 <- ne_countries(type = "countries", country = "Germany", scale = "medium", returnclass = "sf")
map2 <- rnaturalearth::ne_states("Germany", returnclass = "sf")
p1 <- ggplot(map1) + geom_sf()
p2 <- ggplot(map2) + geom_sf()
p1 + p2
library(geodata)Warning: package 'geodata' was built under R version 4.4.2d <- worldclim_country(country = "Jamaica", var = "tmin", path = tempdir())
terra::plot(mean(d), plg = list(title = "Min. temperature (C)"))
library("chirps")Warning: package 'chirps' was built under R version 4.4.2location <- data.frame(long = 100.523186, lat = 13.736717)
d <- get_chirps(location, dates = c("2020-01-01", "2022-12-31"), server = "ClimateSERV")writing: substituting ENGCRS["Undefined Cartesian SRS with unknown unit"] for missing CRSFetching data from ClimateSERV Getting your request...ggplot(d, aes(x = date, y = chirps)) + geom_line() +
labs(y = "Precipitation (mm)")
library(rnaturalearth)
library(elevatr)Warning: package 'elevatr' was built under R version 4.4.3elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'. Use
of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being
deprecated; however, get_elev_raster continues to return a RasterLayer. This
will be dropped in future versions, so please plan accordingly.library(terra)
map <- ne_countries(type = "countries", country = "Switzerland", scale = "medium", returnclass = "sf")
d <- get_elev_raster(locations = map, z = 9, clip = "locations")Mosaicing & ProjectingClipping DEM to locationsNote: Elevation units are in meters.terra::plot(rast(d), plg = list(title = "Elevation (m)"))
library(osmdata)Warning: package 'osmdata' was built under R version 4.4.3Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyrighthead(available_features())[1] "4wd_only" "abandoned" "abutters" "access" "addr" "addr:*" head(available_tags("amenity"))# A tibble: 6 × 2
Key Value
<chr> <chr>
1 amenity animal_boarding
2 amenity animal_breeding
3 amenity animal_shelter
4 amenity animal_training
5 amenity arts_centre
6 amenity atm placebb <- getbb("Barcelona")
placebb min max
x 2.052498 2.228356
y 41.317035 41.467914For example, we can obtain the hospitals in Barcelona by specifying its bounding box placebb and using add_osm_feature() with key = "amenity" and value = "hospital" as follows.
hospitals <- placebb %>% opq() %>%
add_osm_feature(key = "amenity", value = "hospital") %>%
osmdata_sf()Motorways can be retrieved using key = "highway" and value = "motorway".
motorways <- placebb %>% opq() %>%
add_osm_feature(key = "highway", value = "motorway") %>%
osmdata_sf()library(leaflet)
leaflet() %>% addTiles() %>%
addPolylines(data = motorways$osm_lines, color = "black") %>%
addPolygons(data = hospitals$osm_polygons,
label = hospitals$osm_polygons$name)library(wbstats)
indicators <- wb_search(pattern = "poverty|unemployment")
# print(indicators)d <- wb_data(indicator = "MO.INDEX.HDEV.XQ", start_date = 2011, end_date = 2011)
print(head(d))# A tibble: 6 × 9
iso2c iso3c country date MO.INDEX.HDEV.XQ unit obs_status footnote
<chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr>
1 AO AGO Angola 2011 47.7 <NA> <NA> <NA>
2 BI BDI Burundi 2011 48.5 <NA> <NA> <NA>
3 BJ BEN Benin 2011 53.2 <NA> <NA> <NA>
4 BF BFA Burkina Faso 2011 45.8 <NA> <NA> <NA>
5 BW BWA Botswana 2011 80.3 <NA> <NA> <NA>
6 CF CAF Central African … 2011 32.9 <NA> <NA> <NA>
# ℹ 1 more variable: last_updated <date>library(rnaturalearth)
library(mapview)
map <- ne_countries(continent = "Africa", returnclass = "sf")
map <- dplyr::left_join(map, d, by = c("iso_a3" = "iso3c"))
mapview(map, zcol = "MO.INDEX.HDEV.XQ")library('spocc')Warning: package 'spocc' was built under R version 4.4.3df <- occ(query = "Bradypus variegatus", from = "gbif",
date = c("2000-01-01", "2019-12-31"),
gbifopts = list(country = "CR"),
has_coords = TRUE, limit = 1000)
d <- occ2df(df)library(sf)
d <- st_as_sf(d, coords = c("longitude", "latitude"))
st_crs(d) <- 4326
mapview(d)